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Landauer's formula relates the conductance of a quantum wire or interface to transmission prob- 
abilities. Total transmission probabilities are frequently calculated using Green function techniques 
and an expression first derived by Caroli. Alternatively, partial transmission probabilities can be 
calculated from the scattering wave functions that are obtained by matching the wave functions in 
the scattering region to the Bloch modes of ideal bulk leads. An elegant technique for doing this, 
formulated originally by Ando, is here generalized to any Hamiltonian that can be represented in 
tight-binding form. A more compact expression for the transmission matrix elements is derived and 
it is shown how all the Green function results can be derived from the mode matching technique. 
We illustrate this for a simple model which can be studied analytically, and for an Fe|vacuum|Fc 
tunnel junction which we study using first-principles calculations. 

PACS numbers: 73.63.-b, 73.40.-c, 73.20.-r, 85.35.-p 



I. INTRODUCTION 

Since the discovery of giant magnetoresistance in 
metallic multilayers there has been considerable interest 
in studying electronic transport in layered materials 
At the same time, experimental control of the lateral 
scale has enabled studies of electronic transport in quan- 
tum wires of atomic dimensions j 3 " Because of the small 
dimensions involved, the transport properties of both of 
these systems should be understood on the basis of their 
atomic structure. This perception has generated a large 
effort in recent years to calculate the conductance of mul- 
tilayers and quantum wires from first principles. Several 
different approaches have been formulated which have a 
common basis in the Landauer-Biittiker approach or are 
equivalent to it. In the linear response regime, the con- 
ductance G is expressed as a quantum mechanical scat- 
tering problem* and can be simply related to the total 
transmission probability at the Fermi energy, T(Ep), as 

g = p ^t(e f ). (i) 

The multilayer or quantum wire is generally considered 
as a scattering region of finite size, sandwiched between 
two semi-infinite ballistic wires. Aiming at a materials 
specific description, most current approaches treat the 
electronic structure within the framework of density func- 
tional theory (DFT), 5 ! 6 ! 7 ! 8 ! 9 ! 10 ! 11 ! 12 ! 13 : 14 

Frequently the conductance is calculated using a Green 
function expression first derived by Caroli et oZi 15 i 16 
An alternative technique, suitable for Hamiltonians that 
can be represented in tight-binding form, has been for- 
mulated by Ando. 17 It is based upon directly match- 
ing the wave function in the scattering region to the 
Bloch modes of the leads. The latter technique has 
been applied to conductance calculations at the empiri- 
cal tight-binding levelji 8 - as well as on the first-principles 
DFT leveli 8 i 19 i 20 i 21 The relationship between the mode 



matching 17 and Green function 1516-22 approaches is not 
immediately obvious. Indeed, it was recently stated that 
Ando's approach is incomplete and does not yield the 
correct expression for the conductance^ 



In this paper we demonstrate that the two approaches 
are completely equivalent. In the Green function ap- 
proach, a small imaginary part must be added to or sub- 
tracted from the energy in order to distinguish between 
the retarded and advanced forms. 5 ! 6 ! 7 ! 10 ! 11 ! 12 ! 13 ! 14 : 22 In 
mode matching, scattering wave functions are calculated 
that incorporate the retarded or advanced boundary con- 
ditions directly. This makes it possible to solve the 
scattering problem also at real, instead of complex en- 
ergies. In addition to yielding the total conductance, by 
focussing on wave functions the contribution of each indi- 
vidual scattering channel can be identified. In particular, 
we derive a simple, compact expression for the transmis- 
sion matrix elements, see Eq. <|67|) . 



The paper is organized as follows. In the next section 
the Hamiltonian we will use is introduced. This model 
allows us to study both quantum wires that are finite 
in the directions perpendicular to the wire, and systems 
that are periodic in these directions such as single inter- 
faces, sandwiches and multilayers. We will use the single 
term "quantum wire" to describe both systems. In Sees. 
IIIII and llVl the mode matching and Green function tech- 
niques are summarized. The equivalence of the trans- 
mission matrices obtained using these two techniques is 
demonstrated in Sec. Ivland the Caroli expression for the 
conductance is derived from the mode matching expres- 
sions. In Sec. IVII the two techniques are applied first to a 
simple analytical model, 23 and then to an Fe|vacuum|Fc 
tunnel junction using numerical first-principles calcula- 
tions. The main conclusions are summarized in Sec. IVIII 



II. HAMILTONIAN 
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We set up a tight-binding representation of the Hamil- 
tonian. This is not a severe restriction since a first- 
principles DFT implementation that uses a representa- 
tion on an atomic orbital basis set has the same math- 
ematical structure as a tight-binding modelM^ Alter- 
natively, an implementation that uses a representation 
of the Hamiltonian on a grid in real space, can also be 
mapped onto a tight-binding models We begin by divid- 
ing the system into slices ( "principal layers" ) perpendic- 
ular to the wire direction^ The thickness of these slices 
is chosen such that there is only an interaction between 
neighboring slices. Labelling each slice with an index i, 
the Schrodinger equation of the quantum wire becomes 

-H M _ lCl _! + (El - H M ) c< - H M+1 c l+1 = 0, (2) 

for i = — oo, . . . , oo. Assuming that each slice contains N 
sites and/or orbitals, is a vector of dimension N con- 
taining the wave function coefficients on all sites and/or 
orbitals of slice i. The N x N matrices H^j and Hj^ii 
consist of on-slice and hopping matrix elements of the 
Hamiltonian, respectively. I is the N x N identity ma- 
trix. A schematic representation of the structure of the 
Hamiltonian is given in Fig. ^ 

Eq. (J2J is valid both for quantum wires that are finite 
in the directions parallel to the slices, and for layered 
systems that are periodic in these directions. In the lat- 
ter case, translations in the transverse direction can be 
described in terms of a Bloch wave vector in the two- 
dimensional Brillouin zone, ky , which is a good quantum 
number and the system becomes effectively one dimen- 
sional. Explicit expressions for the Hamilton matrix el- 
ements depend upon the particular localized orbital ba- 
sis or real space grid representation usedi2& Since details 
of the tight-binding muffin tin orbital scheme used in 
Refs. [8.19.20 2jJ are given in Ref. |25| and of the real- 
space high-order finite difference method can be found 
in Ref. l2rl they will not be discussed further here. 

The system is divided into three parts, with i = 
— oo, . . . , corresponding to the left lead (L), i = 1, . . . , S 
to the scattering region (S) and i = S + 1, . . . , oo to 
the right lead (R). The leads are assumed to be ideal 
wires characterized by a periodic potential. It is then 
sufficient to identify a slice with a translational period 
along the wire. By construction, the Hamilton matrix is 



the same for each slice of the leads, i.e. = H 



H 



L/Rj 



i_i = B L/R and H M+ i = Bl , R for the left/right 
leads. Fig. n summar i z es our model of a quantum wire. 



III. MODE MATCHING APPROACH 

Eq. J5J) can be solved in a particularly convenient way 
by a technique which we call mode matching. In this, 
we follow Ando and generalize his approach to a slice 
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FIG. 1: (Color online) Hamilton matrix of a quantum wire 
divided into slices. The left (L) and right (R) leads are ideal 
periodic wires that span the cells i — — oo,...,0 and i — 
S + 1 , . . . , oo , respectively. The scattering region spans cells 



A. Bloch matrices 

The first step consists of finding solutions for the leads, 
for which Eq. can be simplified to 



-Bc,_! + (EI - H) a - B+Ci+i = 0. 



(3) 



These equations hold for i = — oo,...,— 1 and i = 
S + 2, . . . , oo, i.e. the left and right leads. To keep the 
notation as simple as possible, we have omitted the sub- 
scripts L and R, see Fig. ^ Since the leads are peri- 
odic wires, one can make the ansatz that the solutions 
have Bloch symmetry, i.e. c, = Ac;_i, Cj+i = A 2 c.;_i, 
where A is the Bloch factor. Substituting this into Eq. © 
results in a quadratic eigenvalue equation of dimension 
N . The latter can be solved most easily by transforming 
it to an equivalent linear (generalized) eigenvalue prob- 
lem of dimension 2N and solving this by a standard 
algorithm. 29 30 

It can be shown that the equation generally has 2A^ 
solutions, which can be divided into N right-going modes 
and N left-going modest labelled as (+) and (— ) in 
the following. Right-going modes are either evanescent 
waves that are decaying to the right, or waves of constant 
amplitude that are propagating to the right, whereas left- 
going modes are decaying or propagating to the left. We 
denote the eigenvalues by A„(±) where n = 1,...,N, 
the corresponding eigenvectors by u n (±) and write the 
eigenvalue equation as 

-Bu„(±)+(£I - H) A„(±)u„(±)-B+A n (±) 2 u„(±) = 0. 

(4) 

In the following we assume that the vectors u„(±) are 
normalized; note that in general they are not orthogonal. 

One can distinguish right- from left-going evanescent 
modes on the basis of their eigenvalues; right going 
evanescent modes have |A(+)| < 1 and left going evanes- 
cent modes have |A(— )| > 1. Propagating modes have 
|A(±)| = 1, so here one has to determine the Bloch ve- 
locity and use its sign to distinguish right from left prop- 
agation. The Bloch velocities are given by the expression 
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V n (±) 



2a 



Im [A„(±)u„(±)tBt u „(±)] 



(5) 
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where a is the slice thickness. A derivation of this equa- 
tion is given in Appendix A. 

Since the eigenvectors are non-orthogonal, it is conve- 
nient to define dual vectors u„(±) 



U n(±) U m(±) = <5n,m 5 u n(±) u m(±) = <W 



(6) 



Any wave function in the leads can be expressed as a 
linear combination of the lead modes. This can be done 
in a compact form using the two N x N Block matrices 
for right- and left-going solutions 



N 



F(±) = ^A„(±)u„(±)Gt(±). 



(7) 



For any integer i, F l is given by a similar expression but 
with A n replaced by A^ . From the foregoing it is easy to 
show that the Bloch matrices obey the equation 



-BF _i (±) + (EI H) B'F(±) = 



(8) 



A general solution in the leads can be expressed in 
terms of a recursion relation 

Ci = C,( + ) + Ci(-) 

= F^(+)c J (+)+F^(-)c J (-). (9) 

Fixing the coefficients in one slice then sets the boundary 
conditions and determines the solution for the whole lead. 



B. 



Transmission matrix 



The scattering region is defined by i — 1, . . . , S, see 
Fig. ^ Right and left of the scattering regions one has 
the recursion relations for the states in the leads from 
Eq. © 

c_x = F- 1 (+)c (+) + F- 1 (-)c (-) 

= [F^(+) - F^(-)] co(+) + F^RccOLO) 
with Cp = cq(+) + cq(— ) and 



cs+2 = F R (+)c s+1 (+) + F R (-)c s+1 (-), 



(11) 



where the subscripts L and R distinguish between the 
Bloch matrices of the left and right leads. 

The boundary conditions for the scattering problem 
are set up in the usual way. The vector Co(+) is treated 
as the source, i.e. as a fixed incoming wave from the left 
lead. There is no incoming wave from the right lead, so 
we set cs+i(— ) = 0. 

Having set the boundary conditions, Eqs. I|1U|) and 
can be used to rewrite Eq. @ in the region not covered 
by Eq. Q, i.e. for i = 0, . . . , S + 1. This describes the 
wave function in the scattering region and the matching 
to the solutions in the leads. Eq. J2J in this region is 
rewritten as 

-H^-iCi-i + ( EI - h m) c * - H^ +1 c l+1 = Q 4 c (+), 

(12) 



with a modified Hamilton matrix defined so H' ; ■ = Hij 
for the elements {i,j = 0, 1}, {i = 1, . . . , S ; j = i, i ± 1} 
and {i, j = 5 + 1, S}, but with 



H' 



Hl + BlF^R, 



h s+i.s+i - H R + B R F R (+). 



(13) 



H' 



and FF 



S+l.S+2 



0, so the modified scattering 



region is decoupled from the leads. On the left-hand side 
of Eq. I|12fl . we have a source term with 



Qo = B L [F^f+J-F^H] 



(14) 



and Qi = for i = 1, . . . , S + 1. Eq. |Q defines a set 
of (S + 2) x N linear equations. Because the Hamilton 
matrix is block tridiagonal, each block being of dimension 
N, this set of equations can be solved efficiently using 
a block Gaussian elimination schemed The total wave 
function c, can then be obtained by back substitution. 

The transmission is obtained from the wave function 
in the right lead cg + i(+). In particular, choosing the 
incoming wave as one of the modes of the left lead, i.e. 
cq(+) = UL, m (+), generalized transmission matrix ele- 
ments T„ jrn are defined by expanding cs + i(+) in modes 
of the right lead 



N 



cs+i(+) 



H ur >'' 

n=l 



( + K 



(15) 



By letting cq(+) run over all possible incoming modes of 
the left lead UL, m (+); m = 1, . . . , N, a full transmission 
matrix is obtained. 

Matrix elements can be defined for all modes, prop- 
agating and evanescent, but of course only matrix ele- 
ments where n, m denote propagating modes contribute 
to the real physical transmission. These modes can be 
selected by making use of their eigenvalues; see the dis- 
cussion following Eq. (0} . The physical transmission ma- 
trix elements are then found by normalizing with respect 
to the current 



' ^R,n(+)«L 

VL, m (+)an 



(16) 



where fL,m(+) and WR jrl (+) are the Bloch velocities in 
the direction of the wire for the right-propagating modes 
to and 77 in the left and right leads, respectively, see 
Eq. JSJ; cil and or, are the slice thicknesses of left and 
right leads^ The total transmission probability is given 
by 



T(E) 



(+) 



(17) 



and the conductance is given by Eq. Q evaluated at 
E = Ep. 
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C. Green function matrix 

Solving the set of linear equations Eq. (|12(l directly 
leads to the conductance. However, to facilitate a con- 
nection to the Green function approach discussed in 
the next section, we can formulate the solution in a 
slightly different way. A finite Green function matrix 
G- ■ (z), i, j = 0, . . . , S + 1 can be defined by 

— H- i _ 1 G-_ lj + (zl - H- Vt ) G- • - H- i+1 G' i+1 j — IS it j, 

(18) 

with z complex. Note, however, that the matrices Hq 
and H^ +1 s+1 are non-Hcrmitian and Gy(£) is also 
uniquely defined for real energies. The Green function 
matrix allows the solution of Eq. (|12[1 to be written as 



= G' 4 



ii0 (B)Qoco(+). 



(19) 



As before, the transmission can be extracted at i = S- 
and comparison with Eq. (fT5|) gives 



U R,„( + )G' S+ l,o(£)QoUL, m (- 



(20) 



which can be used in Eq. (jlfijl . The Green function ma- 
trix block Gg i j (E) can be found by solving Eq. i|18|) 
using a recursive algorithm that resembles a Gaussian 
elimination scheme i 2 ^ 33 



IV. GREEN FUNCTION APPROACH 

An apparently quite different route to the transmis- 
sion matrix starts by defining an infinite Green function 
matrix Gi,j(z) for i, j = — oo, . . . , oo with respect to the 



original Hamiltonian of Eq. J2J ■ 



H, < |G, + (zl — II, , 



Gj,j — H^j+iG 



(21) 

Choosing z = \im v ^o(E ± irf) defines as usual the re- 
tarded/advanced Green function matrix. We shall use 
Gij(E) to denote the retarded Green function matrix 
and G"j (E) to denote the advanced Green function ma- 
trix. 



Here gh(z) and g R (z) are the surface Green functions of 
the semi- infinite left and right leads, respectively, which 
can be calculated using an iterative technique: denoting 
G[™! (z) as the solution of an equation similar to Eq. pijl. 
but with Hi,j = for {i > n V j > n}, one can easily 
derive the right-going recursion relation 

zl — H n _|_i !n _|_i — H n+ i : „G^ ri (z)H„ ! „ + i 

G&SU(*) = L (23) 

For an ideal wire with i, j = —oo, . . . , n, G[™k(z) = gL(-z) 
should be independent of n resulting in the following 
equation for the surface Green function, 



zI-H L -B LgL (z)B[l g L (z)=I 



(24) 



Several iterative algorithms have been formulated for 
solving this non-linear matrix equation^2^*^& A simi- 
lar reasoning based upon a left-going recursion for the 
right lead results in an equation for the surface Green 
function gR,(z) of the right lead 



zI-Hf 



B^g R (z)B R 



gR.0) 



(25) 



Again, setting z = E + irj in Eqs. I|24|) and (|25|l defines 
the usual retarded surface Green functions gL/R (-£-)■ Al- 
though we are mainly interested in the physical limit 
lim^Oj m practice a finite value of r\ is often used in 
order to make the iterative algorithms stable. 
The quantities 

£ L (£) = B LgL (£)B[ ; E R (£) = B R g R (£)B R , (26) 

which appear in Eqs. (|22|I - H25|I . are called the self- 
energies of the left and right leads, respectively. 16 Once 
these are obtained, the finite Hamilton matrix of Eq. i|22|) 
is constructed and the retarded Green function matrix 
Gij (E) can be found using a recursive algorithm^ Us- 
ing the lead modes the transmission matrix elements can 
then be calculated, as will be shown in the next section, 
Sec. lIVBl Alternatively, the total transmission probabil- 
ity can be expressed in a form that does not require the 
lead modes explicitly, which is discussed in Sec. IV1 



B. 



Transmission matrix 



A. Partitioning 

Eq. (|21|) is most conveniently solved by applying a par- 
titioning technique. 16 - 34 It is straightforward to show that 
the finite part Gj j (z) defined for i, j = 0, . . . , S+l can be 
derived from a closed set of equations, similar to Eq. (|21|l . 
but with Hij replaced by H" ■ where H" 3 = H^- for the 
elements = 0, 1}, {i — 1, . . . ,S; j — i,i ± 1} and 
= S +1, 5}, but with 



H£ (z) = H L +B LgL (z)B[, 



H 



S+1.S- 



.!(«) = H R + B R g R (z)B R . 



(22) 



The transmission matrix can be obtained from the 
Green function matrix of Eq. I|21l) . To do this, we 
adapt a Fisher-Lee type of approach to our tight-binding 
formulation. 37 Assuming that the unperturbed reference 
wave function is the Bloch mode UL, m (+) that comes in 
from the left lead, the Lippmann-Schwinger equation 38 
in tight-binding form is 



c <-./ v /> i| i- 



m . k 



( + ) 



fi(+) + ]Tg m v jV£ f£(+) 



-)■ (27) 
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Here UL,m,i(+) is the reference wave function in slice i. It 
obeys Bloch symmetry and u Ljm (+) = u Lim! o(+) ls the 
Bloch mode at the origin, see Sec. inTH The matrix V,- k 
represents the perturbation with respect to the ideal left 
lead. 

Eq. (|27J) can be simplified using the Dyson equation, 
which in tight-binding form reads 

Gi, = G^p + Gj,3'Vj,fcGfc°o 

0,k 



fL(+) + E g ^' v ^ f l(- 

3,k 



G& (28) 



using Eq. Comparing Eqs. (O and (PBJl one finds 

the simple expression 



Ci = G ifi (E) G$l(E) u L , m (+) 



(29) 



From the definition of the generalized transmission ma- 
trix elements, cf. Eq. (|15() . one then obtains the expres- 
sion 

Tn,m = UL„ (+)G S+lfi (E) \G$(E)] 



U L ,m(+). (30) 



To find T n ^ m one needs to calculate only the Green func- 
tion matrix blocks Gs+i,o(-B) of the full system and 

Gq°q(.E) of the ideal left lead. The physical transmission 
matrix elements and the total transmission probability 
can then be obtained from Eqs. (|16J) and (|17(l . 



V. MODE MATCHING VERSUS GREEN 
FUNCTIONS 

The two seemingly different formalisms introduced in 
Sees. IIIII and lTvl are in fact closely related. In this section 
we will show how all Green function results can be ob- 
tained from the mode matching approach. We begin by 
expressing the Green function matrices of ideal wires in 
terms of the Bloch matrices, F(±). These expressions are 
then used to prove that the transmission matrix elements 
obtained from the mode matching and Green function ap- 
proaches, cf. Eqs. (|2(J[1 and l]3Up. are in fact identical. Af- 
ter that, we show that the transmission matrix elements 
are independent of the exact positions within the leads 
that are used to match the leads to the scattering region, 
apart from a trivial phase factor. Then we derive from 
the mode matching expression for the total transmission 
probability the Green function expression known as the 
Caroli expression^ Finally, a more compact expression 
for the transmission matrix elements is derived. 



A. Green functions of ideal wires in terms of Bloch 
matrices 

We begin by deriving an expression for the retarded 
Green function matrix G^j of an ideal infinite wire in 



terms of its eigenmodes. The columns of such a Green 
function obey the equation 



[E+I-H) Gj?-BT G ^. = I^, 



(31) 



where E + = E + irj. For i ^ j the solution is simi- 
lar to that of the wave functions, see Eq. J2J. In ad- 
dition, the retarded Green function should consist only 
of propagating waves that move outwards from the <5- 
source and/or evanescent states that decay away from 
the source^ From Eq. (0, we have the column recur- 
sion relations 



G^j(E) = F*-'(-)G^(E),i<j, 



G\°](E) = F i ->(+)G<fj(E),i>j 



(o), 



(32) 



The diagonal block GfJ (E) can now be obtained by com- 
bining Eqs. l|31fl and (|32[1 . which gives for i = j 



[-BF- X (-) + E+I - H B t F(+)] Gf] = I. (33) 



Eliminating E + I — H using Eq. (JHJ) then yields 



or the equivalent 



G<°>(25) "=Bt[F(-)-F(+)] 



B 



F- 1 (+)-F- 1 (-) 



3,3 



(34) 



(35) 



Eqs. (|32|l and l|34|l represent the full expression for the 
Green function Gf} of an infinite ideal wire in terms 
of the Bloch matrices F(±) and thus in terms of the 
eigenmodes. Note that we can set E + = E since, in 
terms of the modes, the retarded Green function matrix 
is uniquely defined for real energies. 

The advanced Green function matrix g| q Q (-E) can be 
found from a similar procedure. It should consist of 
propagating waves that move towards the source and/or 
evanescent states that grow towards the source. One can 
construct two new Bloch matrices F a (±), which are sim- 
ilar to those defined in Eq. (JJJ). In F a (+) one collects 
the modes that are decaying to the right (growing to the 
left) and modes that are propagating to the left. F a (— ) 
then contains modes that grow or propagate to the right. 



F Q (±) 



N 



n=l 



A£(±K(±K f (±), with 



(36) 



K(T), <(±) = u„(=f) propagating 
A„(±), u^(±) = u„(±) evanescent. 



Using these definitions, expressions for the advanced 
Green function matrix are obtained from Eqs. l-il'l) and 
(GU) by replacing F(±) with F a (±). 

From the general relation between retarded and ad- 
vanced Green functions, Gjj = (G^)^, the following 
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row recursion relations can be deduced for the retarded 
Green function 



i < J, 



G?J(E) = G<$(E) [F a H-)Y~ l , i> j. (37) 



The retarded Green function G\ s } (E) of a semi-infinite 
wire extending from i = — oo, . . . , can be obtained using 
a similar technique. Instead of Eq. I|32|) . we get 



G, ( ; ) (£),PHg(E), 1 <0. 



(38) 



where g(E) — Gq S q(-E) is the surface Green function. 
Using this in Eq. I|31l) gives for i = and j — and for 
i = — 1 and j = 0, respectively, 

[-BF- 1 (-)+£;+I-H]g = I, 

[-BP- X (-) + B+I - H] F- 1 (-)g = Btg. (39) 

Note that the term is absent in the first equation 
since we are dealing with a semi-infinite wire. These two 
equations can be easily solved to find an expression for 
the surface Green function 



g(i?)=F- 1 (-)(Bt)- 1 . 



(40) 



Eqs. (|38|l and l|4U|) represent the Green function of a semi- 
infinite ideal wire extending from i = — oo, . . . , 0. In 
a similar fashion, one gets for the Green function of a 
semi-infinite ideal wire extending from i = 0, . . . , oo 



Gl s >{E)=F t+1 {+)B-\i>0. 



(41) 



Analogously to Eqs. (|38l) - (|41|l . one can define the ad- 
vanced Green function matrix G[ S Q a (E) in terms of the 

Bloch matrices F a (±). Moreover, since [g a ]^ = g, we 
have the following relation between the Bloch matrices 



B t F(±) = F at (±)B. 



(42) 



B. Equivalence of mode matching and Green 
function approaches 

The retarded surface Green functions of the left and 
right leads can be derived from Eqs. I|4(J|) and (|41H 



g L ( J B) = F I : 1 (-) 



b; 



; Sk (E)=F k (+)B^. (43) 

The retarded self-energies of Eq. (|26[) are then given by 

E L (£) - BlF^H ; Sr(£) = B^F R (+). (44) 
Comparing Eqs. I|13|) and (|22H then establishes 

G' i<j (E) = G iij (E). (45) 



In other words, the two Green functions discussed in Sees. 
11111 and HVI are identical. 

By comparing Eqs. (|14|> and (|34|l one has 



Qn 



(46) 



In conclusion, the two expressions for the generalized 
transmission matrix elements, Eqs. I|20[l and (|30|l . are 
identical. 



C. Invariance of transmission probability 

Apart from trivial phase factors, the transmission ma- 
trix elements should not depend on where in the ideal 
lead the wave function matching is carried out. In a re- 
cent paper it was stated that Ando's expression for t n _ m , 
Eqs. I|2UI) and (|lfc>|) . lacks this invariance property and is 
therefore incomplete. 22 One can however prove directly 
from Eq. (|20|l that the transmission matrix elements do 
have the required invariance property^ The proof be- 
comes easier if the equivalence of Eqs. (|20fl and 13U|I . 
established above, is used. 

The scattering region runs from slices to S + 1 if 
we include the boundaries with the left and right leads. 
This means that the Green function matrix Gi_j with 
indices i,j outside this region obeys the equations for the 
ideal leads. From the column and row recursion relations, 
Eqs. lj3"2")> and (|3*7|l . one derives 



Gs+i+ij{E) - F R (+)G S , +1 JE) 



(47) 



for j < 0, i > 0. In a similar way, one derives for the 
Green function matrix of the left lead 



G?](E) - Fi(+)G$(25) FJT(-) 



(0), 



r a t/ 



(48) 



for j < 0. 

We now artificially extend the scattering region by in- 
cluding slices from the left and right leads and let it run 
from j < to S + 1 + i with i > 0. Eq. gives for the 
transmission matrix element 



G$UE) u L>m (+) 



T 'n,m - U R^n( + ) G S+l+i,j( E ) 

Using g7J) and gSJ) then gives 

<,m = 5L(+) F k(+)G s+1 , (i?) 



(49) 



F L J (+)u L , m (+) 



= A R ,„(+)A-M+K, 



(50) 



with a real. The third line follows from applying Eq. J7J . 
The last line follows from the fact that we are only in- 
terested in propagating states and for propagating states 
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A | = 1. Using this result in l|16|) and l|17fl proves the 
invariance of the total transmission probability with re- 
spect to moving the boundaries between leads and scat- 
tering region into the leads. 



D. The Caroli expression 

The total transmission probability is given by Eq. (|17fl . 
where the sum has to be over propagating states only. We 
can extend the summation to include all N states (propa- 
gating and evanescent) by defining an N x N transmission 
matrix 



t = V^)rV L 5 (|), 



(51) 



where r is the matrix whose elements are given by 
Eq. (Hi. Vr(+) is defined as the singular, diagonal ma- 
trix that has the velocities i>r jTI times the constant ft/<ZR 
on the diagonal for the right-propagating states and ze- 
ros for evanescent states. We call it the velocity ma- 
trix. Likewise a pseudo-inverse velocity matrix Vl(+) 
can be defined, which has I/wl,™ x a-h/h on the diagonal 
for left-propagating states and all other matrix elements 
are zero. These velocity matrices project onto the space 
of the propagating states so that the transmission ma- 
trix has only non-zero values between propagating states. 
(|17|l can then be expressed in the familiar form 



T = Tr[th] 

= Tr pV R (+)rV L (+) 



(52) 



Using the above definition of the velocity matrix and 
Ando's expressions for the transmission matrix elements, 
we will show how (|52ll can be rewritten as 



T = Tr [r R GT L G°] , 



(53) 



where G r ,G a are short-hand notations for Gs+i,o{E) 
and G% S+1 (E), respectively. The matrices r L / R are de- 
fined as 



L/R 



£ L /R - S L/R 



(54) 



Eq. l|53|l is known as the Caroli expression^ and it is 
often used to calculate transmission probabilities ^LSSi It 
is equivalent to the Kubo-Greenwood expression for the 
linear response regime 

The first step is to construct N x N matrices U(±), 
the columns of which are the eigenmodes u„(±), and 
diagonal matrices A(±), the elements of which are the 
eigenvalues A„(±) 



U(±) = (ui(±) u 2 (±) ... u N (±) 



(55) 



From Eqs. © and © it is then easy to show that the dual 

vectors u n (±) form the columns of the matrix [U(±) _1 ] 
and that the F(±) matrices obey the equation 



F(±)U(±)=U(±)A(±). 



(57) 



In a similar way matrices U a (±) and A a (±) can be con- 
structed, see Eq. (PEJ- 

Using these definitions, one can generalize the r-matrix 
of Eq. JU to 



r = U R 1 (+)G f 'Q U L (+). 



(58) 



Note that r-matrix elements are defined not only between 
propagating states, but also between evanescent states. 
However, as we remarked above already, only propagat- 
ing states contribute to the physical transmission. 

The second step is to express the velocity matrices in 
terms of the T-matrices. To do this we use an expression 
for the velocity matrix, 

V(±) = i [U t (±)B t U(±)A(±) - A t (±)U t (±)BU(±)] , 

(59) 

which can be shown (see Appendix A for a proof) to be 
equivalent to the definition introduced in the first para- 
graph of this section. Using (|57|l , this can be rewritten 
for the right lead as 

V R (+) = iU R (+) [b r f r (+) - F R (+)B R 1 U r (+) 



= *u R (+) 



Sr - St 



Ur(+) 



u R (+)r R u R (- 



(60) 



The second line follows from l|44|l . A similar relation 
between the r-matrix and the velocity matrix for the 
left lead can be shown to exist, 



v L (+) = u^(-)r L u£(-), 



(61) 



by using (|42[l and an equivalent expression for the ve- 
locity matrix, Eq. (|A7|) . Eqs. $§U$ and JHU imply that 
the T-matrices project onto the space spanned by the 
propagating states. 

The third step is to introduce a matrix P that explic- 
itly projects onto the propagating states of the left lead 



p = u L (+)i p [u£Hr 



N„ 



(62) 



where the I p -matrix contains 1 on the N p diagonal ele- 
ments that correspond to propagating states, and at all 
other positions. Given this projector matrix, it is possi- 
ble to prove that 



A(±)r 



A„(±)£„ 



(56) 



Q P = ir L . 



(63) 
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The proof is given in Appendix B. Using this property 
one has 

Q U L (+)V L (+) = Q PU L (-)V L (+) 

= *r L u£(-)v L (+) = i [uf (-)l 1 V L (+)V L (+) 



P (3 Ji L P R P (3 



T t 



Ul T (-) 



(64) 



-00, 



* = D, 1, 2, 



. X: 



Substituting JSSJl, JSDl an d <EH i nto HU leads directly 
to the Caroli expression, Eq. (JHSJ- 



E. Transmission matrix: a compact expression 

Using the results of the previous section it is possible 
to derive a more compact expression for the transmission 
matrix. Combining (|51|) . I|58|) and i)62[l one has 



FIG. 2: The parameters of a one-band nearest neighbor tight- 
binding model for a single impurity in a onedimensional chain. 



The roots A(±) can be given a more familiar form. For 
|oj| < 1 we define a wave number k by 



cos(fca) = w, 



(70) 



where a is the lattice parameter. From Eqs. (|69[) . H7U|) 
one then obtains 



t = V R (+)U R 1 (+)G''QoPU£(-)V L (- 



(65) 



A(±) 



(71) 



By following the same steps as in Eq. 164(1 and using 

~ i i 

VlV£ = this can be simplified to 

t = *v|(+)U R 1 (+)G''[U^(- ) ]-i V J(+). (66) 

Writing out the transmission matrix elements gives the 
compact expression 



which describes propagating states. For \oj\ > 1 one de- 
fines k by 



cosh(fta) 



(72) 



V aRflL 

This is the tight-binding equivalent of the Fisher-Lee 
expression relating transmission and Green function 
matrices i^i 



One obtains A(±) = exp(=F«a) if lj > 1 and A(±) — 
— exp(=FKa) if a; < —1; both cases describe evanescent 
states. 

Since the scattering region consists of a single im- 
purity, S — 1, Eqs. lfT2"|> - lfl4T give three linear equa- 
tions with three unknowns describing the scattering prob- 
lem. In a one-channel model one has F(±) = A(±) and 
A(±) _1 = A(^). There is only one possible incoming 
wave, so co(+) = 1. The linear equations then become 
in matrix form 

















-( 


VI. EXAMPLES 









/3[A(-)-A(+)] 





(73) 



A. A simple analytical model 

We consider a system consisting of a single impurity in 
a one-dimensional chain and treat this within a one-band 
nearest neighbor tight-binding model. The parameters of 
this model are given in Fig. [21 This model can be solved 
analytically, 23 so it can serve as a simple test to illustrate 
the equivalence of the different approaches. 

We will first solve the problem using the mode match- 
ing approach of Sec. lIIII It is convenient to define a scaled 
energy by 



with 



E — e — /3A(+) -f3 L 
A= I -/3 L E-d 

-/3r E-e- /3X(- 

' (74) 

Solving this set of equations and using Eqs. lj7U|) and i|7T|) 

we obtain the compact expression 



C2 



2ika 



-if sin(fca) 



d + (1 — b) cos(fca) — i6sin(fca) ' 
defining the dimensionless parameters 



E 



2/3 



(68) 



Pl + Pl d . e-Q ± /?iA 



2/3 2 



20 



, / 



ft 2 



(75) 



(76) 



The model involves only one channel and with u m (±) = 1 
Eq. J3J reduces to 



Applying Eqs. (|15|) - (|17f) then yields for the total trans- 
mission probability 



-/3+(£-e)A(±)-/3A(±) 2 = 0, 



(69) 



T(E) = \c 2 \ 



(77) 
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Using Eqs. I|70|) and l|75[l it is easy to show that this trans- 
mission probability is identical to Eq. (15) of Ref. |U 
which was obtained using a different technique. 

It is instructive to solve the same problem using the 
Green function approach of Sec. IIV Al First one has 
to find the surface Green functions of the leads from 
Eqs. (|24(l and l|25|l . which for the current model become 



[E — e — [3 2 g(E)] g(E) = 1, 



(78) 



where E is a real energy. This equation has the solutions 



9 ± (E) 



(79) 



for both leads. The Green function matrix in the scat- 
tering region can then be found from Eqs. l|21|l and Q22JI. 
which can be combined in the 3x3 matrix equation 



AG(-E) 



(80) 



where Gij(E), i, j = 0, . . . , 2 are the matrix elements of 
G(E) and A is given by Eq. lj7"4T) . Inverting A yields the 
matrix element 



f 



2/3 d + (1 — b) cos(ka) — ib sin(fca) 1 



(81) 



with the parameters b, d and / defined by Eq. (|76|) . The 
relevant Green function matrix element for the ideal lead 
is found from Eq. (|34J) 



2/3sin(fca)' 



(82) 



Using these results in Eq. (|30(l one observes that the ex- 
pression for the (one channel) transmission matrix ele- 
ment becomes identical to Eq. (|75|l . 

Finally one can calculate the transmission probability 
from the Caroli expression given in Sec. lVDl cf. Eq. (|53|l . 
Using Eqs. J2BJ), |2} and GHJ one obtains 



= — 2/3sin(fca), 



(83) 



for left and right leads. Using Eqs. (JH2>, (E3 an d 
Gq 2 = (G2,o)* m Eq. (|53l) then yields an expression for 
the transmission probability that is identical to Eq. (|77|l . 
It illustrates the equivalence of the different approaches 
for calculating the transmission in this simple model. 

In addition to providing a channel for propagating 
states, an impurity can also give rise to localized states, 
whose energy is outside the energy band of the chain, cf. 
Eq. (|72|) . Such a state does not contribute to the phys- 
ical transmission, but the transmission amplitude has a 
pole at an energy that corresponds to a localized stated 
Within the mode matching approach this corresponds to 
an energy at which c$+i becomes infinite. For the present 
model the energies of localized states can be obtained by 



setting k — in and setting the denominator to zero in 
Eq. J7SJ). This leads to the equation 

(w + d) (u + sgn{uj)y/oj 2 - l) -6 = 0; \u\ > 1, (84) 

the roots of which give the energies of the localized states. 
Again these results are equivalent to the results obtained 
using the approach of Ref. |U 

Within the Green function approach the energies of the 
localized states are given by the poles of the Green func- 
tion matrix. Via Eq. H81fl this again leads to Eq. H84fl . 
Alternatively, since the Green function matrix is the in- 
verse of the A- matrix, cf. Eq. (jHUJ, its poles are given 
by the roots of det(A) = 0. This equation is equiva- 
lent to Eq. 184[) . as is easily shown by setting A(+) = 
±exp(— Ka) in the A-matrix. 



B. Fe|vacuum|Fe tunnel junction 

As an example of a more complex system, we consider 
an Fe|vacuum|Fe tunnel junction where the electronic 
structure is treated using the local density approxima- 
tion of DFTii£ The calculations are based upon a tight- 
binding muffin tin orbital (TB-MTO) atomic spheres ap- 
proximation (ASA) implementation£ii2i22i2ii2£ of the for- 
malism described in Sect. IIIII 

The first step in the calculation is the self-consistent 
determination of the electronic structure of the tun- 
nel junction using the layer Green function approach of 
Ref. |3f|. The Fe leads are oriented in the (001) direc- 
tion and the atoms at the Fe(001) surfaces are kept at 
their unrelaxed bulk positions. For the bcc structure 
and TB-MTOs, 4 i a principal layer in the (001) direction 
contains two monolayers of Fe with a thickness of 2.866 
A. The vacuum region is modeled by a number of such 
slices, of the same thickness, filled with "empty" spheres 
of the same size and packing as the Fe atomic spheres. 
The atomic sphere potentials of the vacuum region and 
four monolayers (two principal layers) of Fe on cither 
side of the vacuum are calculated self-consistently while 
the potentials of more distant layers are kept at their 
bulk values. These potentials then form the input to a 
transmission calculation based on mode matching&22i2i 
Further technical details can be found in Ref. |2!| 

A useful quantity to extract from the self-consistent 
layer calculation is the layer density of states (LDOS) 
Pi(E). It is related to the retarded Green function matrix 
defined in Eq. (|2*T)l by 



p. t (E + i7j) 



-7r _1 ImTr[G M (£; + z7?)] (85) 



where the trace refers to the usual Im angular momentum 
indices characterizing MTOs. For reasons of numerical 
stability the retarded Green function matrix is calculated 
by adding a finite imaginary part to the energy; we have 
used -q = 0.0025 Ry. In the mode matching approach, 
the LDOS can be directly expressed in terms of the wave 
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FIG. 3: (Color online) Layer density of states (LDOS) at the 
Fe(001) surface layer. Top panel: majority spin. Top curve 
(red, dashed): as calculated using the layer Green function 
methodi-- Bottom curve (blue, solid): as calculated using the 
mode matching approach. Both calculations use r\ = 0.0025 
Ry. For clarity, the top curve has been displaced by 10 units 
along the y-axis. Bottom panel: minority spin, the LDOS is 
shown with a negative sign. 



functions or, alternatively, in terms of the Green function 
matrix of Eq. ifTHj l 



Pi {E) = -TT^ImTrtGUS)]. 



(86) 



This Green function matrix can be calculated for a real 
energy, but in order to make a comparison to the results 
obtained with Eq. Q85|l. we add an imaginary part, r\ — 
0.0025 Ry. 

In Fig. |21 we compare the LDOS obtained using 
(|85|l and (|86() for the topmost Fe monolayer of (001) 
Fe|vacuum|Fe where the vacuum layer was so thick (four 
principal layers, corresponding to a thickness of 11.466 
A) that the LDOS corresponds closely to that of a free 
Fe(001) surface. The two curves, displaced vertically 
for clarity, are indistinguishable, illustrating the equiv- 
alence of the two Green functions defined in (|18fl and 
(|2"TJl . Moreover, the LDOS are in essential agreement 
with results found previously for a Fe(001) surfaced 

If one resolves the LDOS into contributions from dif- 
ferent parts of the surface Brillouin zone, then the con- 
tribution at r (k|| = 0) exhibits sharp peaks near the 
Fermi level. These are associated with characteristic 
surface states found on (001) surfaces of bec transition 
metals, 45 These surface states are mainly derived from 
d 3z 2_ r 2 orbitals on the surface atoms projecting into the 
vacuum, and they have Ai symmetry4S^ They become 
most clearly visible if one filters out the contribution to 
the LDOS at T of the states with Ai symmetry. 

This contribution, calculated using the mode match- 
ing approach, is shown in Fig. 0] as a function of the 
distance from the surface. In the majority spin LDOS, a 
sharp peak is found at Ep — 1.8 eV, which is very promi- 
nent in the surface Fe layer, cf. Fig. Eld) and decays 




J 








-2 -4 -2 • 

EC/*?) 



FIG. 4: (Color online) The top panels represent the band 
structures along T-H of bulk Fe for majority and minority 
spins. Panels (a-h) give the Ai contribution to the LDOS 
pi{E) at r of the Fe| vacuum interface as a function of the 
layer position, (a) LDOS of an Fe layer at an infinite distance 
from the interface; the double peak structure below -4 eV 
and single peak (of a double peak structure) above -1 eV 
correspond to the Ai bulk bands, (b) LDOS of the 6th Fe 
layer (the surface layer is layer 1). (c) LDOS of the 3rd Fe 
layer; the sharp peak in the bulk band gap corresponds to 
a surface state, (d) LDOS of the the surface Fe layer; the 
surface state peak has maximum amplitude, (e-h) Same for 
minority spin. 
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FIG. 5: Calculated kii -resolved transmission in (a) the major- 
ity and (b) the minority spin channels through a Fe-vacuum- 
Fe tunnel junction. The vacuum region has a thickness of 
5.733 A and the Fe electrodes have a parallel magnetiza- 
tion. Only the central part of the Brillouin zone is shown, 
i.e. — 7r/3a < k x ,k y < ir/3a; the transmission in the other 
parts is close to zero. 

rapidly into the bulk, cf. Figs. Elc)-(a). This evanes- 
cent state clearly represents a state that is localized at 
the surface. In the minority spin LDOS, a sharp peak 
with similar properties is found very close to Ep. It 
corresponds to the surface state that is observed experi- 
mentally by STM4^ The positions of both these surface 
state peaks in the majority and minority LDOS are in 
excellent agreement with those obtained using a Green 
function (KKR) approachi 4 ^ 4 ^ Such surface states are 
good examples of evanescent states. Clearly, they are 
described properly within the mode matching approach 
which was, however, disputed in Ref. 22. 

Using the mode matching approach, we next calculate 
the transmission through an Fe|vacuum|Fe tunnel junc- 
tion in which the vacuum region is only 5.733 A thick 
(corresponding to two principal layers) and the mag- 
netizations are parallel. Fig. [5] shows the calculated 



ky-resolved transmission for majority- and minority-spin 
channels. The majority conductance has a maximum at 
T and it decreases smoothly going away from T. This be- 
havior is quite general for scattering from a potential bar- 
rier. At fixed energy, a wave travelling perpendicular to 
the vacuum barrier penetrates furthest into the vacuum 
and therefore has a maximum tunneling probability^ 

In contrast, the minority conductance is dominated by 
sharp spikes of very high intensity, close to 0.13 FX. This 
behavior has been analyzed in terms of interface reso- 
nant states that extend relatively far into the vacuum. 
States originating from the two Fe| vacuum interfaces cou- 
ple through the thin vacuum region, which enhances the 
transmission considerably. 4 — These transmission spikes 
or "hot spots" are observed quite generally in calcu- 
lations for perfect transition metal | insulator | transition 
metal tunnel junctions42i§£ A calculation using Eq. i|30|) 
gives the same transmission within numerical accuracy. 
It again illustrates the correct treatment of evanescent 
states within the mode matching approach. We have also 
tested expression (|67|) numerically and found that it too 
gives the same transmission within numerical accuracy. 

Finally, the invariance property of the transmission dis- 
cussed in Sec. IV CI has been tested numerically by mov- 
ing the interfaces between the leads and the scattering 
region, where the matching is carried out, into the leads. 
In this way more and more slices of lead are treated as 
scattering region, see Fig. ^ the transmission should be 
invariant under this operation. Adding up to ten princi- 
pal layers (20 monolayers) of bulk Fe to each side of the 
scattering region changes the transmission negligibly, by 
less than one part in 10 8 . 



VII. SUMMARY AND CONCLUSIONS 

We have demonstrated the equivalence of the mode 
matching and the Green function approaches to calcu- 
lating the conductance of quantum wires and interfaces. 
In the mode matching technique, the scattering problem 
is solved by matching the wave function in the scattering 
region to the Bloch modes in the leads. The technique 
is formulated for tight-binding Hamiltonians and covers 
all representations that can be expressed in tight-binding 
form, including first-principles implementations using lo- 
calized orbital basis sets or a real space grid. 

Alternatively, the scattering information can be ex- 
tracted from the Green function, which is calculated by 
partitioning the system into a scattering region and leads. 
We demonstrate that the Green function technique can 
be reformulated in terms of mode matching. In addi- 
tion we prove that the mode matching expression for the 
transmission matrix does not depend on where in the 
leads the scattering region is matched to the ideal leads, 
which was called into question in Ref. I22I 

Calculating the full transmission matrix allows us to 
study individual scattering amplitudes from and to ev- 
ery possible mode. We have derived a compact expres- 



12 



sion for the transmission matrix elements, Eq. I|67[) . Only 
propagating modes enter this expression, since evanes- 
cent modes do not contribute to the transmission directly. 
The evanescent modes are important, however, for set- 
ting up the correct matching conditions at the bound- 
aries between the scattering region and the leads. Alter- 
natively, the total transmission probability can be cal- 
culated from the Caroli expression, Eq. (|53|l . Formally 
this expression sums over all possible modes, i.e. propa- 
gating and evanescent. However, by deriving the Caroli 
expression from the mode matching approach we show 
explicitly that only the propagating modes give a non- 
zero contribution. 

The mode matching approach and its equivalence to 
the Green function approach were illustrated for a sim- 
ple analytical model, as well as by numerical calculations 
on an Fe|vacuum|Fe tunnel junction. In the latter we 
treat the electronic structure within DFT, using a TB- 
MTO-ASA basis set. We demonstrate that the layer den- 
sities of states that follow from the mode matching and 
Green function approaches are numerically indistinguish- 
able. We identify the Fe(OOl) surface state in the density 
of states and establish its contribution to the transmis- 
sion through the tunnel junction. Finally, the invariance 
property discussed above is demonstrated numerically on 
the Fe|vacuum|Fe tunnel junction. 



APPENDIX A: VELOCITY MATRIX 



In this appendix we will prove that 



V(±) 



-v m (±)5„ 



m,n — ty m\- 1 -/ u m,n; 

a 



for the expression introduced for the velocity matri- 
ces in Sec. IV Dl Eq. Q59p. Here, a is the translation 
period along the wire and, for right- respectively left- 
propagating states, u n (±) is the Bloch velocity in the 
direction of the quantum wire. For evanescent states 
v n {±) = 0. 

For ease of notation we drop the index ± in the follow- 
ing. From Eq. Q and its complex conjugate one has 

-u m Bu n + X n u m (EI-H)u n -(X n ) 2 u m B^u n = 0, 
-u^Bt^ + \* m ul (EI - H) u„ - (A^) 2 ut,Bu„ = 0. 



Multiplying the first equation by X m , the second by A„ 
and subtracting the two gives 

[A„uJ„B t u n - A^u^BuJ (1 - X* m X n ) = 

-iV m ,„ (1 - A;„A„) = 0, (A3) 

according to Eq. So if A* n A„ ^ 1 then V m ,„ = 0. 

The velocity matrices contain, by construction, either 
right- or left-going modes. For right-going evanescent 
modes one has |A| < 1, for left-going evanescent modes 
A | > 1, so A* n A„ ^ 1 if n and/or m denotes an evanescent 



mode; thus V m ^ n = 0. For propagating states, |A m | = 1 
and one can write A m = exp(ifc m a). It is obvious that if 
k m i= k n , then A* n A„ ^ 1 and again V TOjIl = 0. 

This argument does not hold for the diagonal matrix 
elements of propagating states, and also not for degener- 
ate propagating states, i.e. if X m = X n . In the latter case 
we take the derivative d/dE of the first line of Eq. 1A2|) . 
The coefficients of the terms in du n /dE and du^/dE 
vanish because u„ and uj^ satisfy the eigenvalue equa- 
tion Eq. |0J and its complex conjugate, respectively, and 
A* = 1/A„ for propagating states. Collecting the remain- 
ing terms gives 

|X u L Btu n - A*u^Bu„] + A„u|„u„ = o 



dE 



(A4) 



where we emphasize that for m ^ n this only holds for 
degenerate propagating states. Without loss of gener- 
ality, degenerate states can be chosen orthogonal, i.e. 
u^Un = for m ^ n. Since dX n /dE ^ 0, it then follows 
that V m „ = if m ^ n. 

The diagonal matrix elements m = n for propagating 
states are also given by Eq. (|A4|) . Setting n — m and 
using the fact that the states are normalized, i.e. u}jU„ = 
1 gives the expression 



V -iX ™- 

v n.n — t/x n , , 

dX n 



(A5) 



Since for propagating states we can write A„ = exp(ika) 
and the Bloch velocity is defined as v n = h~ 1 dE /dk, this 
(Al) then proves both Eq. IjAlfl and Eq. 0. 

It is straightforward to derive an alternative expres- 
sion for the velocity matrices in terms of the advanced 
matrices 

V(±) = i [U Qt (T)B t U Q (T)A a (T) - 

A»t(T)U"f(T)BU a (T)] . (A6) 

Multiplying from the left by [A a ^] 1 and from the right 
by [A a ] _1 this expression is seen to be equivalent to 

V(±) = i [{A at (T)} _1 U Qt (T)B t U Q (T) - 
(A2) U a t( T )BU a ( T ){A a ( T )}- 1 ] . (A7) 



APPENDIX B: PROJECTOR MATRIX 

In order to prove Eq. (|63[) . we start from Eq. (|14|l and 
rewrite it using Eq. I|42l) as 



Qo = BlFl - Fjt(-) 



B 



(Bl) 



In the following we will drop the subscript L for ease of 
notation. Multiplying (JBlfl on the left with the identity 
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operator I = X) u «(+) u n(+) an d on the right with I 

E<(-K, f (-) y ields 



N 



(B2) 



An(+) 



<t(-)Bu„(+) 



a*t(-) 



<t(_) B t Un (+) 



where use has been made of Eqs. (JBJ, @ and 136|l . 

By similar arguments as those leading to Eq. (|A3() it 
follows that [. ..] = 0, unless \? n {-) = 1/A*(+). From 
Eq. (|36| it follows that this is true if n = m and n denotes 
a propagating mode. In addition, for every (+) evanes- 
cent mode n with eigenvalue A„(+), there is one (— ) 
evanescent mode with eigenvalue A ra (— ) = l/A*(+)i2i 
Again it follows from Eq. H36fl that only the m = n 
terms are non-zero for evanescent modes in Eq. (IB2|I . 
This equation then simplifies to 



Q = E<(-)St( + ) 



(B3) 



A„(+) 

N 



ut (+)Bu„(+) - A„(+X(+)Btu„(+) 



n=N p + l 
1 



A„(+) 



ul RBu„(+) - A„(+)ut (-)Bt Un (+) 



where the first summation is over the N p propagating 
modes and the second summation is over the N — N p 
evanescent modes. 

The projector matrix P is defined by Eq. (|62[1 . In 
the product QoP only the first summation in Eq. I|B3|) 
survives. Moreover, since [. . .] = iV n ,n according to 
Eq. (|A3|) . we have 

N p 

Q P=igu»(-Kt(-)V n , n (+). (B4) 

71=1 

Together with Eq. ||HTJ this then proves Eq. (JHSJ- 
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